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ABSTRACT 

“This report describes the development and implementation of a computer 
algorithm, based on the shallow-fluid equations of oceanography, for calculating 
+ the wind field above complex terrain. The algorithm contains a two-dimensional 
shallow-fluid model in the form of a fully-documented computer program compatible 
with a UNIVAC 1108 machine. To guide the selection of initialization procedures 
and the optimum finite-differencing scheme applicable to numerical solutions of the 
algorithm, a detailed study was made of the analytical solution of the shallow-fluid 
equations for one-dimensional flows over an isolated ridge. The selected optimum 
finite-differencing procedure is a Lax-Wendroff scheme using nine grid points and 


two time levels in combination with a nine-point low-pass filter. a 


~ 


A comprehensive computational program, using an isolated symmetrical 
mountain, was carried out to provide guidelines as to the nature of two-dimensional 
solutions of the shallow-fluid equations for the wide variety of initial conditions 
encountered in the atmosphere. The model was initialized by impulsively accelerat- 
ing the fluid to a constant velocity everywhere and, after a while, the flow near the 
mountain approaches a steady state. The results showed that the flow patterns ' 
could be divided into four major categories: Subcritical without hydraulic jumps; 
supercritical without upstream waves; critical with hydraulic jumps and wind- 
direction reversals; critical with hydraulic jumps but without wind-direction rever- 
sals. For the subcritical and supercritical flows, initialization procedures do not 
appear to pose a problem. For the critical flows, care must be taken in the selec- 
tion of initialization procedures. 


Comparisons of calculated wind field patterns with recent detailed observa- . 
tions of wind circulations above complex terrain show excellent qualitative agree- 
ment in the limited cases available for analysis. Additionally, the computer ue 
algorithm for the two-dimensional model, when applied to the terrain at White ——~_ 


>> Sands Missile Range, gave results that were consistent with limited observations 
available for two example situations. 


The computer program containing the two-dimensional] shallow-fluid model 
is written in FORTRAN V language and is fully documented in the four appendices 
to the report. The documentation includes user's instructions, a complete program 


listing, detailed flow diagrams, and a completely worked example problem. 
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USER INSTRUCTIONS FOR THE ASL/WSMR WIND FIELD TERRAIN 
ADJUSTMENT MODEL COMPUTER PROGRAM 


The ASL/WSMR Wind Field Terrain Adjustment Program is designed to 
adjust the wind speed and direction field to the terrain heights at each coordinate 


intersection over a reference grid system. 


The computer program is written in the FORTRAN V language and is 
designed for use on a UNIVAC 1108 computer. The program requires 31390, 0 words 
of core storage including system and FORTRAN library routines. The computer 
program has the optional capability of printing the adjusted wind speed and direction 
fields and surface layer heights at specified time intervals. In addition, the wind 
fields and layer height information may be written on magnetic tape for subsequent 
graphical processing. The program uses FORTRAN Logical Tape 5 (Card Reader) 
and Logical Tape 6 (Printer) for standard input/output. The program also uses 
Logical Tape 1 as a file for optional output. The computer program construction 
is described ta Section A. 1 and user instructions are given in Sections A.2, A.3 and 
A.4. An example problem is given in Section A. 6 and a program listing is shown in 
Appendix B. : 


A.1 COMPUTER PROGRAM CONSTRUCT 


The Wind Field Terrain Adjustment Model Program consists of six aub- 
routines including the main driver program MODEL. Program MODEL determines 
the number of problem cases to be executed in a single computer run and determines 
the core allocation of variable program arrays for each problem. The computer 
program uses object time-dimensioning techniques to accommodate different grid 
system sizea. The grid system is assumed to be the first quadrant of a Cartesian 
coordinate syatem with the positive y-axis oriented north and the positive x-axis 
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oriented east. The total size of the grid system matrix is limited by the 
equation: 
N 2 7*IDIM*JDIM+2*IDIM+2* JDIM 


The parameter IDIM is the maximum size of the x dimension and JDIM is the 
maximum size of the y dimension of the grid system matrix. The parameter N is 
the DIMENSION'd size of a large array Q in MODEL containing all program variable 
arrays that depend on the grid system matrix. The parameter N is currently set to 
18000, but it can be set to any value in order to accommodate the program to com- 
puter core storage limitations. 


The first subroutine, Subroutine JACK, is the main calculation routine. 
This subroutine inputs the model data, determines the program options desired, 
sets the initial wind field conditions and then enters a time loop through the shallow- 
water model equations discussed in Section 4 of the main body of the report. The 
computer program uses a second-order central space and forward time differencing 
scheme to evaluate Equations (4-14) through (4-21). The wind speed and direction or 
the u and v components of the wind speed, as well as the surface layer heighta, are 
printed and/or written to tape at selected points in the time loop. When the time _ 
loop is completed, the program stops or loops to the next problem in sequence, 


The remaining subroutines are utility support programs for Subroutine 
JACK, Subroutine MISC determines the area of uniform grid spacing within the 
grid system and returns the indices of this area for tape output. Subroutine UVDIR 
converts the u and v components of the wind speed to wind speed and direction for 
printing. Subroutine MOUTNR Inputs the terrain data into a specified area of the 
grid system. The program fille any remaining area of the grid with height values 
calculated from the minimum height and the height at the edge of the input area as 
a function of distance. Subroutine OUTPT is the tape output routine and uses the 
UNIVAC NTRAN routines. 
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A logic diagram of the computer program is given 'n Section A.§ and 
detailed flow charts of the computer program are given in Appendix D. 


A.2 PROGRAM INPUT PARAMETERS 


The data input parameters required for the computer program are listed 


in Table A-1. The information categories in the table are defined as follows: 


CARD GROUP 
SEQUENCE 
NUMBER - 


NAMELIST - 


FORTRAN - 


MODEL - 


UNITS - 
LIMITS - 


VALUE Oe 


ARRAY SIZE 7 
CARD COLUMNS - 


FORMAT - 


Order of input of the three card groups. 


Name of the FORTRAN NAMELIST list to which the input 
parameters belong (Card Group 2 only). 

FORTRAN symbolic notation defining the program input. 
Mathematical notation correspondirg to the FORTRAN 
notation. 

Dimensional units of the input parameters. 


Numerical limits on input parameters. 


Default value substituted if the present value is zero (Card 
Group 2 only). 


Maximum number of core locations reserved for the parameter. 
Data card punch field. 


FORTRAN input format. 
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A.3 DATA INPUT METHOD 


This computer program uses formatted as well as namelist input state- 


/ ments. The parameters using a formatted read statement are self explained in 


Table A-1. The namelist input data must be in a specific form in order to be read 


YY 
using a NAMELIST list. The first character in each card to be read must be blank. 


The first card in the namelist list contains the namelist name preceded by the 


character $. The last card in each namelist list contains $END to terminate the 


list. The form of the remaining data items in the list may be: 


a. Vartable Name = Constant - The variable name may be a sub- 


scripted array name or a single variable name. Subscripts must be integer con- 


stants. The conetant may be integer or real. es 


b. Array Name = Set of Constante (Separated by Commas) - The array 


name is not subscripted. The set of constants consists of constants of the type 
The number of constants must be less than or equal to the array 


real or integer. 
Successive occurrences of the same constant can be represented in the form 


size. 
k* constant. 


The sequence of the input data parameters within the list is not significant. 
A more detailed explanation of the FORTRAN NAMELIST can be found in any Fortran 
The input parameters within the namelist that have default values 
Parameters 


Language Manual. 
in Table A-1 are initialized to zero prior to input of the first case. 


that are not used or have default values need not appear in the namelist list. When 


multiple cases are stacked, all parameters retain their values from the previous 


case and are changed only by input. 


A.4 EXPLANATION OF PROGRAM INPUTS 
The program input parameters are arranged into four card groups. 
A.4.1 Card Group Number 1 


This card group contains information specifying the number of 
cases to be executed and the size of program storage arrays. 


ND - Number of times to loop through the entire program reading a new set 
of values for NP, IDIM and JDIM on each loop. 


NP - Number of times to loop through the main calculation routine using 
constant values of IDIM and JDIM and reading Card Group 2‘(also Card 
Group 3 if requested) on each loop. 


Maximum number of grid coordinates in the x direction for all cases 
within the NP loop. 


IDIM 


JDIM - Maximum number of grid coordinates in the y direction for all cases 
within the NP loop. 


(IDIM and JDIM are used to determine the size of program storage 
arrays at execution time. See Table A-1 for the limits on these 


parameters. ) 
A. 4.2 Card Group Number 2 
This data card group contains most of the program input data. 


The data parameters in this group are read using the Fortran Namelist QLST1 to 
simplify program input and reduce the number of input cards for similar cases 


stacked in sequence. Table A-1 gives default values for applicable parameters in 


this card group. Also, all parameters in this card group retain the value set by the 


previous case unless changed by input. 


LL - Number of grid coordinates in the x direction for the present case. 
JL - Number of grid coordinates in the y direction for the present case. 
xX - Array containing the coordinates of the x-axis of the grid system in 


ascending order. 


Y - Array containing the coordinates of the y-axis of the grid system in 


ascending order. See Table A-1 for default values for x and y. 


ISKIP Program option control flag. All ISKIP options are assumed initially 


zero, 


a. If ISKIP(1) is set to 1, the calculated wind field and layer height 


values are printed at each time given in the array PRINT below. 


b. ‘If ISKIP(2) is set to 1, the wind field and layer height arrays are 
output to tape at each time given in the array PRINT. If ISKIP(1) 
and ISKIP(2) are both zero, then ISKIP(2) is set to 1. 


c. If ISKIP(3) is set to 1, the terrain height matrix is output to tape 
(ISKIP(2) must equal 1 also). 


d. If ISKIP(4) is set to zero, the program assumes ULT and VLT 
below are input as the initial u and v components of the wind speed 
and the adjusted u and v components are printed in the wind field 
output print file. 


If ISKIP(4) is set to 1, the program assumes ULT and VLT are 
input as the initial u and v components of the wind speed and the 
adjusted wind speed and direction are printed in the wind field 
output print file, 


ABLK 


PRINT 


ULT 


VLT 


If ISKIP(4) is set to 2, ULT and VLT are assumed to be the 
initial wind speed and direction and the adjusted u and v com- 
ponents of the wind speed are printed in the wind field output print 
file. 


If ISKIP(4) ia set to 3, ULT and VLT are assumed to be the initial 
wind speed and direction and the adjusted wind speed and direction 
are printed in the wind field output print file. 


e. If ISKIP(5) is set to 0, the terrain height matrix is input in Card 
Group Number 3. 


If ISKIP(5) is set to 1, the terrain height matrix used in the last 


executed case is used in the present case. 


. 


Minimum layer depth. If the layer depth is found to be less than ABLK 
in any iteration, the layer depth is then set to ABLK. A value of 50 


meters is assumed if zero is input. 


Array containing the iterative time step value in minutes at which the 
wind field and layer height values are to be printed and/or output to 
tape. Values are arranged in ascending order. 


The initial u component of the wind speed if ISKIP(4) is set to 0 or 1. 
The initial mean wind speed if ISKIP(4) is set to 2 or 3. 


The initial v component of the wind speed if ISKIP(4) is set to 0 or 1. 
The initial direction if ISKIP(4) ts set to 2 or 3. 


The initial surface layer height. 
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DTLMDA - Stability factor used in calculating the time step DT. This value should 
be as close to 1 as possible. The program assumes a general value of 
0.95 if zero is input, Program instability is indicated by a program 
message IFLAG ERROR, When this occurs, either the value of ! 
DTLMDA or NCNT must be reduced. 


v p 0 
Gl - Reduced gravity factor, where Gl = s(} - #) re c( - #] 
0 1 
= the acceleration of gravity 
0, = the potential temperature at the top of the surface layer 
i @ on the potential temperature at the bottom of the surface layer 


P4»P) = respective layer densities 


The program assumes a value of 0.1 for Gl if zero is input. 


ISMOTH - The number of time steps between the application of a nine-point filter 
to the adjusted values of the wind field. A value of 10 is assumed if zero 
is input, If no smoothing is desired, sct ISMOTH to a large value, say, 
10, 000. 


The number of time steps between the recalculation of the time step 


time increment DT. A value of 2 is assumed if zero is input. 


The Fortran logical output unit for tape output. The program uses 
NTRAN L/0 with unit 0 an illegal unit. The program assumes unit 1 if 
zero is input. If more than one reel of output is executed, a second reel 


must be assigned to IUNIT + 1. See Section A.7 for the output tape format. 
A.4.3 Card Group No, 3 


This card group contains the terrain heights and the starting and 
ending indices that define the area of storage of the terrain heights within the grid 


system. The terrain heights outside of the input areaaare calculated from the 


minimum height and the height on the edge of the input terrain. The starting and 

ending indices of the x and y axes are input first, followed by the terrain heights. 

IST - Starting index of the terrain height matrix in the x direction. 

IND - Ending index of the terrain height matrix in the x direction. 

JST - Starting index of the terrain height matrix in the y direction. 

JND - Ending index of the terrain height matrix in the y direction. 

HG - The matrix of terrain heights input by the following Fortran statements: 
DO 10 J= JST, JND 


10 READ (5, 2000) (HG(I, J), I=IST, IND) 
2000 FORMAT (14X, 10F6, 1) 


A.5 SUBROUTINE LINKAGE FOR THE COMPUTER PROGRAM 


The logical linkage for the computer program subroutines is shown in 
Figure A-1. Each connector represents a communication link between the 


| subroutines. 
A.6 EXAMPLE COMPUTER PROGRAM 


This section explains the input data example shown in Figure A-2. The 
example consists of a problem with WSMR terrain elevation data. The example 
data shown here is a computer listing of the data. The data are stored in a program 


file as an element of the file. 


The first card image shown is a system (@MAP, }) card starting in column 
one. This card is used to direct the system to link the program subroutines and 
form an absolute program deck. The second card image @XQT directs the system 
| to load and execute the program using the data that follow. 


A-12 


FIGURE A-1. Diagram of linkage between subroutines of the computer program. 


*eyep ynduy atdurexa aq; jo 3uj38}] Jaindulog = °Z-y JZUNDLI 


E°nOs fe VOZTs OTZT I SHETSOTSZS HGS IV’ Hen G OPC IT Inc TL aret fer Sornne 
@°e7ut an, sS°re@ 
@°LOwTS* UZOTS°UZHTL® BTOST COL ey HGten Ls CPU UBSEH°LeCZ9 COTS fat esneeo 
O°GSOTO°SOST a “WLZTO*SOSTS HT cl2°EtZic *Zeet7?°otssrt eres oret foc 26nree 
S*alcto°shste OtZta° BLETS*SSCeI*OSSTE wtntICcUSETO°Sceto’Ssest flon 150998 
e°er7et fon eterna 
O° ROUT OCSTS° GOTZT °90616° 06079 ° HE TET* HOTZ9° CETL OCTZE OL te 900 6%enge 
O°alUteOBels * tOCre*Ssazteezere sezte’szzte’zeztz°orett caret 00 erreee — 
TOTEIO°LOTTs “SOSTH°LELTOOSE SSL ute °ERETT VCE IY ecets OTCt Cor £%O688 
2° teet 600 grenes 
4° OSOTS°OTOTYSZOTE “eSOTS*Z2hO2S °ScEts OSETL* TTUZL UGHT9°2L02 Oav ste0e0 
O°GOSTT* SOmTy OSZTL “Gaste HEete wEZiC*eeztzolZiz 9212 °otet 800 qrORnee 
TETSETCTSTE°STITE°osOTS*ZL0Zn* Lents ONCTE SELIG *SzET6 e7Et Geo evonee 
o°sest Ano seeeee 
O°L9LTO°L9L Te LELI9°ORGTS° 7220726 Q6ETL*OGLI0* Te91e°Szeteezet 000 1867900 
O°90LT9° SCTE QOCTO LOSTE *SZZte *uretS LETTE Seziz*oyerrt eter 600 eveese 
ViSTEDTSESRELT ISTO SOLIS" SERIE’ SOETE*OOE TO UPELE SCETE SRE! seo ecocee 
0°azs 100 
O° OZUT a LCLTU°ZS9T6° 9OLTO°SZETO“SOGTL *CCOTS UZE16°SnOTe OnSt Gao 
Q°GMST ET ZONTyOTET O° veZTE *SEZto* NOZIT*ZONThH HESIE °ezzte otal 6e0 
TET TU" LORI OSE TONES TO*EETET “LISTS COTS" SEL TT: Teste ceet 008 
S°Cee 


6°90L Te SHOTS SESTE°DOLTS°LOLTS “GENTE £2516 °LTSIO HZSIS C6et 
COLTS SHOT E THETO°SSZTO*GSZSTZ* OOTTO W2STL GUZIC *SZZTZ° OTST 
Lert eee otererene re enter Sealy eels Sects OSE TE Tee 
TEES NESTE COSTS “GENTS SNC Te OSnTS ESNTT “O9NTG’ OGeTO ESET 
O°eSE TO OTOTE FONT A°SETIOOSTIZ COTTI ecw, “GestS 4czte* vest 
SHRSUS ERE IE PEE TA TOUS tD SORTS  SOLNE SONtE: LNEIS Bette Oeet 
O° SECT LOLAIT° LEST9° BASTA°OTE BE ASE TO SSHT eH’ LONTL SOSTO°RZEt 
C°QO9TOSOSTO°RESTOSEETOSOTTL "GOTTL “GOZTS 9OTTS “Cazte°eest 
; OP SPETL PIREELO TSO ENGI 3 OESES:OOLTE SNH TT NETS: OOE14 O01 
T°nee 
T*THETOSIRTESONTI OTE TS "SFE TS SOCTe*OZnTh LOeTOezstt erst 
V*mGSTE COSTA SERTO"SECTE *9BTTZ* OUTTA ERZTO°ONTIO CHZtS eect 
°SZSTS*GTSIS*LETTT°ZONTS SSEZE*SOLTT°Anats Conte -otzte gett 
oe 7X 3) 
O° RENTS CSOTS HOC TS COETS*ZENIC ately aunt CONTS Conte azgt 
O°GESTE°SSOTS TIBTO°SSETTSECTEGtCleeczty SeztOCuzteerzt 
S*otcto° wOSTs GeZTS*LOETS “CC TZO° SNOT TeCleSSZte Lette entt 
s 9 9 9 

On2s 

*TSLINNE 

*SEuLowST 

*O0RZEi le 

asa 

75s) 

*O0¢ ‘OSTSINIUed 

O°T*TTTssTywSI 

Tass 

Te TH 

t etenge 

atye zeorce 

3¢a%r* = 000 renAnee 

Coletezn/ser-re0 13 

wito W77°1598 


A-~14 


PUES 


avd e22090 3400 W200w ANSWLSAC OV avnuass OI3Ts Only MwSe/asy 


eel eee aa EE <a 


(penupuod) “2-V AuNols 


SMLAYV MOTS a LELTE NOEL POabe TUMET Hee TATE PE oo ay : 7i-nno 
ere th coSty SLU 29 Cem eg te EL be OULTE eh ae  e? itrane 
atrenn 

a TILTH EMU be ISATL AT Ck Tue hE hTERT FO seh eC h esbed 6% renn 
C8QycEV* HOTEsAEHSH SLATE Tye LLO ME ECE TON TEER, Le vitane 
S°AMC TI GNZTS* TLE HEEL CSE’ (ce OZH LL FHS OSG Lh ber TE GHG! Lu-rerg 
reset . at. cAng 

UQUGTE OGwTo LTUZS VE LeL Cres MdCe heh ELF LEG 8 yOne a GLTOAA 
S*AmvsG MLOTY*ECTSL a CaTy Cd spa SG be Pye TS MEL SLL in Teel : a lreng 
“PASS TS LES a SAID * Veg lure2yis yMeer Vers uGHse oo, THe yrgd eo reno 
ys CPyt zutron 

w hGulsOBUTC*ZHUZI EO IPS* Te ce tute es gy Moe Loh LE . torano 
S*A2LUSyOVET LS LILTH®ungle* Thale MHEG HeererueZ le HetL  26eT Oft*no 
S°OMCTE*SSSAa*SSVSTU nLaty*Sr loe*egslo Ogre eLGcd “Tirzh Erat 66"FN0 
2°C%ul 96°00 

L°VWETZS°OCET SS ZHUZST OT 2S ES ecu9GZER°62571 Switch LET Hale L6nang 
S°GSCCS*HLo ta ISETHULSTE LIne HME le *Qn2 ts °anZin one*Gietel : 96°r00 
EPASNSLE* SSSla SST ME ?STSsCoctyoGOLTV G2ETT HILL NULTS 2OGT sénrno 
S°Leut #6000 

C*OSKTR ObBTo LO0ZT*L9T EO CEL eh YEheh OAH 7H ES uct oO HAG COE? eererno 
S*GTa 2S *SMUZO* SHITE HIGTO* Foe TE HZELS CMZTL CGA at teeter! 200700 
CP AlclS°ylelu LOST]°y Lele says le 292 sy LAL TE 292 1G Uk TH ZeYyT tereng 
RPT o6fan0 

bb HCuUlS°O2e Ty 666TR OC ZCL* 19 ZO* LONE GANZ HON H HENCE 1 9GE 6e"ro0 
Y° m6 2S °SS2Z0°S99TOIcoMTL Oe To SZEIS°CO2T2 HhZ ih oct] fuel g8nnno 
C°OTCTUAT2T ur L0STO 2027S *Fbnlb HERI e PATH OMHTO"ECUTT MAST Lenroo 
& oT 9¥"ren 

B°ObwI9ONLTe* TT0ZS°SS729* SLuco OTZcU CHZz%* UC heb* CONZF *N9GZ sennno 
H°MENSYELTZe OLITLo9nIT Toe sy ULC T weZteomZs7] rel terel *¥onoo 
T°OTSAU LOST HOSE *cETIS*LESIn GILG *GHOTE CLLI9 EL TG 22ST eeeroo 
Covel 2anrano 

&°OSUTE*GOLT. “OSOTE*2E0TL OC TAL T9Tcu'T I Hey o2S22"2He2e “T7y92 T¥NNn00 
S°2LuZU°LILTUSGEST2°wOMTO*CES Fu ZEZSTE* HESTS"ONZIZ*E1cTULPel o8eroe 
GHG LIL °SOT Tu TOTS SOT TW GSI Y LOLI LOLTE Y9L EL OC YTC HET 62006 
H°2o4T e2oneo 

B°OSule°ECETS°OZ6TZ*SHucd*Snten* Heeze *“cAnZn* obnZl TILE ° VON! Ldrnge 
O°GSCS9° SLOSS LILT9® iLeTe Oa Te°us[Te*Bu2th Leu le°*oi2Tt etet 9d°NNO 
S°LOLTO* HET To OOZSTL* SOTTO H2nTe G2 h LZCLTY 22219 EUSTC’ tent sernoe 
6°*SagTt edtroe 

H°2EeTL°OS6Tc “OZETT*COT2S *Fyuze* TOEen  GAnZe N9IGCE Ti e2° 2992 eéinnoo 
S*nOuca*eetes OL9TS* Sze le *agztuuSZlo°anzte *zzz1z°o1ztz “ote? zearne 
O°NOCTA* HOTTA OTST EE Teteazylo*GOLIO Zs cTO°SOSIE EL STL Cnt tzor0o0 
zZ°teot oceroe 

Mo mOeTL UTOZe ° TTOSL * TTGSY*CE1C9* SLOTS *O9TESS*OTLEE * Peyet "GOuZ 6 2nr0o 
SO mlw2z Tosh LELTL°9telo*nscleGhZl9°Cn2reeezer2°stety eter €9°700 
B°MGLT6°OOZTL SET TH LEZTO*Ghylt *“GI9OTE*SmOTT*Z2ST0*Cynt2°unnt aoerne 
: ¢°uzgt 990"°e8 

wrR2aTL* LTOSL* TVOST* OC uSO*CE IST Hote esses. *uIGe2°Cuces "LTeS sonnee 
f° NBwZ0° 9BZSY°SL0ZE* B79 OH LY" EHSth HEZte o@lele*yiztoanel ST | 
U°MOLTO° LETT SOTZTO*TL2TO°OOLT TOC USL *CFSTN T2STL AL HTH OMe cvnge 
4 °SFyT zonnne 

B°KSwTL TTOZe* TEETI°TLOST ENT Su *SSZeH*GAN7Z4 OAICS CULTS 29S t99900 
S°ONL29 ECT Su 6GOTZ GhzTh TecIn Tezte* crete UTcle suse ated : 0%9°900 
S°NOTH*OOSTe*OTSTS*HdeTUrs9L lo SHOTL *2NGTHLPHIL OLETO CT osneng 
H°ecat S006 

B°OBETS* 2HG 22° USET9I* L020 *AGT2y*SSZZU* LHEPE LOGS Tue ene 40°00 
S°LEe2L* TUO%e"UOSTTS*ce[|tersecic*otz]peun7z2a*yTZcrerary sce $3°ANO 


€£.u9u 3400 W3U0,. LNAASNCOV LrCUSL Vata Unt MySay 


fad 


RAN I Ek RING ee pe eee edn: SOT he ale at Ra etiy eMbie  aest 


(penuyjuog) =°*Z-Vv ZHNDII 


terol rer scrsro & 
S°ES ben UENZa° 06SZS°SG97TUS Lm susLHesarEs 20° HET e9* Cr 120°C TET fee ecteno & i 
U°L9LTE* GOL 9E9TH BAG TU ZS IS *COmth Line LPSIZ OCG To Gagt on" %Tr99 |S f 
S*2huZ9* L69Te “GESTS USHID*AZHIY* ZEMTL® Tem re L866 9GTE “Cent fan sctraog 4 
O° age ana sc Teno t 
UMOsehLLESu°LOHS9°SLOSL* TIUZL OSG 1° 92070 LET CL UGETH EHOT een ct tag 
S*SSOTS°*LNGTC°HNIDTO*SOSTT*LZS Iu OQMIE "CANTO UIMTH LuMte GayT enw ectego ‘ 
V°GHYTS® TOGT YE SLSTO “COMTI *2Onl es GEMI9*SenTs THM lO OGeTE Ther fer IE Fang i 
G*udet Ll 4 OC °° 90 | 
W°6Oul0°OLOS4°ZL029 °UZOTL PGols *SEHIY* Ss UZO°EEI EO EL TES *YME2 fev 6zTtra0 l 
R°O6n22 ° TKESyOZOTL USETY*SOSTUTEGHIT ZOMTT 20H Ls ELETE HTEt ent @7t"00 
L*TWHIS’CoMld SCTZ0°ugEl9*Zenst Song *Senswe *2tmye Oyt1o* wont fer 4zteng 
: € "erat rer e2Trng 
O° whuth® nOgTu°OZ8TE°9SETS* Tels" OZGI9 CE T79° CE TLS VIE TH OhHe acer szteno 
S*HU529° Hh Cou "EGRTA LITT TESTA LEMIG *2EHTL HPS IL UEC TE Zon ene ezten0 5 
T*SOHIS* ewe LILTL °USETIAIHTy OSHI9ZSEHTN HAUL Tutte “Enekt nes e200 
o*onat te. z2tene 
H°ONLTD ONL Te°9SETH 4ELTOEZETI “CnIZe* CZs Es§led *Sceze*gayz bad @ 4 t2t*e0 
EC WISM*GEnTe° TluZewTETO*AOL ETSONTIL7GTS°LHE LT TANCE zeqt fo ecTeng 
Smtr le Z6LTS°CONTO*SESTE CES TU TISIL 2OS TC lew T OGnTN* yonT “cr 6tt900 
9° 2eut fer etrene % 
MOVIE £92 Tu GUBTZ utETe* 29 TY CCTey Qmlen Seer eut:ze* 1 ty2 Aor Ltrang 
L°MLACE® LE9SU°L9L16°902 TE°Gayly QaHle T/CTN ESS IL Sects ant nel 9treno 
M°AOSTU THE Ty “GUSTS HIGTH SIG le “CEQTS*LOSTA GUS LE LEGTH PGT bal Co ctreao : 
DO HrPgd ent ettene 
VLMa lEGOLTEOZETHESOZE JBL CY RETZY"VLOAT ENT ZE°% CEC 1Z4 he rm ctreno 


oe seve edaguyy 3av0 W360. INPWISWL Kavurse, UTS 2 UT OH egy 


The first data card shows a1 in column two. This value (ND) is the number 


of times the entire program is to be executed reading a new set of dimension para- 
meters on each pass. The second data card sets NP to 1 and IDIM and JDIM both 
equal to 41. The parameter NP represents the number of times the program is to 
be executed using IDIM and JDIM set to 41. All following cards beginning with 
$QLST1 are repeated NP times. The third data card begins the namelist QLSTI. 
The parameter ISKIP indicates: 


(1) The adjusted wind field and layer heights are printed at each time 
given in the array PRINT 


(2) The adjusted wind field and layer heights are output to tape 


(3) The terrain heights are output to tape 


(4) ULT and VLT are input as vector components and mean wind 


speed and direction are printed in the output listing 


(5) The terrain height data will be input following the namelist 
QLST1. 


The remaining parameters in QLST1 indicate that output is to take place at 180 and 
300 seconds iteration time and the u and v components of the wind speed are set to 
7 meters per second. Also, the initial surface layer is 2800 meters high and filter- 
ing is to be done on each fifth time step and the output unit is Unit 1. Several para- 
meters from QLST1 do not appear in the data because default values are provided. 


The parameters and values are: 


LL=41; JL = 41; X = UTM coordinates of terrain (see Table A-1); 


Y = UTM coordinates of terrain (see Table A-1); ABLK = 50 meters; 


DTLMDA = 0. 95; G1 = 0.1; NCNT = 2. 


The namelist is then ended with a $END card. The next card gives the starting and 
ending indices of the x and y axes, respectively, which define the area of storage 
for the terrain heights that follow. The terrain heights begin in Column 15 of the 
card and are read using the statements shown in Section A. 4. 3. 


The above example problem executed in 302 seconds on the UNIVAC 1108 


computer at the University of Utah. The computer output listing for the sample 
problem is shown in Appendix C. 


A.7 OUTPUT TAPE FORMAT 


The output tape produced by the program is a binary (odd parity) tape. The 
data are recorded in integer binary and floating point binary (real) form. The tape 
(tapes) is in multiple file form where each file represents the output from a single 


case. The records of each file are arranged as follows: 
RECORD 1: 


Word 1 - Number of times the arrays UL (u component), VL (v component) 
and PL (layer height) occur within the file (integer). 


Word 2 - Flag where if set to 1 indicates record 4 contains the terrain 
height data HG. If set to zero, then record 4 contains the first occur- 


rence of UL (integer). 


Word 3 - IDIM or the I dimension of the arrays UL, VL, PL and HG 
(integer). 


Word 4 - JDIM or the J dimension of the arrays UL, VL, PL and HG 
(integer). 


; 
Word 5 - LL or the number of values in the x axis in record 2 (integer). 
Word 6 - JL or the number of values in the y axis in record 3 (integer). 
Word 7 - IST or the starting index on the x axis of uniform grid spacing 


(integer). 


Word 8 - IND or the ending index on the x axis of uniform grid spacing 


(integer). 


Word 9 - JST or the starting index on the y axis of uniform grid spacing 


(integer). 


Word i0 - JND or the ending index on the y axis of uniform grid spacing 


(integer). 


| 
! 


q RECORD 2: LL words consisting of the x axis in ascending order (real). 
q RECORD 3: JL words consisting of the y axis in ascending order (real). 


RECORD 4: If Word 2 of Record 1 is set to 1, then Record 4 contains IDIM*JDIM 
words of terrain height data. The data are arranged as if the following 


i. statement were used to write it: 
((HG(I, J), l=1, IDIM)J=1,JDIM) - HG is a real variable. 


RECORD 5: The u component of the wind speed (real). These data are arranged 


as if they were written by the following statement: 


(UL, (I,J, 2), 1=1, IDIM), J=1, JDIM) 


RECORD 6: The v component of the wind speed (real). 


(VL(I, J, 2), 1=1, IDIM), J=1, JDIM) 


RECORD 7: The height of the surface layer (real). 


(PL(I, J, 2), l=1, IDIM), J=1, JDIM) 


A-19 


Records 5, 6 and 7 are repeated (4, 5 and 6 if word 2 of Record 1 is 0) the number of 
times indicated in Word 1 of Record 1. If more than one reel of output is written, an 
end-of-file is placed at the end of the first reel. This file mark will appear within 
the affected file and will indicate a switching of reels is necessary on input. (The 
program will rarely require more than one reel of 2400 foot tape, even with several 
stacked cases.) The program prints the contents of each output tape file at the end 


of each case loop. 
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APPENDIX B 
COMPUTER PROGRAM LISTING 


Appendix B coutains a complete listing of the computer program. The pro- 
gram is written in the FORTRAN V language and has been run on a UNIVAC 1108 


computer. 
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APPENDIX C 
COMPUTER PROGRAM EXAMPLE PROBLEM OUTPUT LISTING 


The computer program output listing as shown here has several pages of 
output omitted because of the volume of the listing. Only important sections have 


been retained. 


The first page of output gives all of the program initial input data except 
for the x and y axes and terrain heights. Pages C-3 through C-18 (computer listing 
pages 25 to 40) give the initial wind field, layer height and terrain values up to a J 
index of 16. This part of the listing was included by setting ISKIP(1) equal to 1. 
Also, by setting ISKIP(4) = 1, the program has printed the wind speed and direction 


rather than the u and v components. Pages C-19 through C-34 (computer listing 


{ 
] 


pages 66 to 81) show a printout of the 79th time step at 3 hours up to a J index of 
16. The listing was produced by setting ISKIP(1) equal to 1. This part of the liat- 
ing gives the x and y coordinates, the wind speed and direction, the layer height 
and the vorticity, where vorticity is a measure of the rotation of the wind in units 
of ‘eo Pages C-35 through C-50 (computer listing pages 107 through 122) print- 
out of the 131st time step at 5 hours up to a J index of 16. This part of the listing . 
has the same form as that given above for time step 79. This part of the listing 
also represents the final solution as the greatest time value input was 5 hours. 
Page C~61 (computer listing page 148) gives a summary of the contents of the out- 
put tape requested. The summary shows that output from two time steps (79 and 
131) was desired and that both were written to tape. The listing then prints the 
time step and number of model seconds. The listing then shows the maximum index 
values of each output array and then gives the indices of the grid area of uniform 


spacing. Also, the summary shows the terrain height data were written to tape and 
the output tape was unit 1. 
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if APPENDIX D : 
COMPUTER PROGRAM FLOW DIAGRAM 


an 
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4 Appendix D contains detailed flow diagrams of the ASL/WSMR Wind Field 
. boos Terrain Adjustment Program. 
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AO ee le 


tm 


bt PROGR MEDEL 


This pregras reade the odject tine dimension limits end calculates the relative addresses of the program errays 
thet depend on grid size aud pacses the addresses to the wain calculation routine JACK. 


(Input the auabder of tines to 
execute entire progres) 


(input che aumber of cizes to 
execute the program using the 
x and y dimensions IDIN aad JOIN) 


READ NP, IDIN ,JDIN 


-. 


Jil « IDDNOIDIN 
332 © JLLe2 


Sie Ja (Calculate the stertiag locetions 
32 = JZ of the variable program errays 
Jé @ J24012 within the erray Q) 

J5 © Jée3Ll 

36 = JStIDIN 

37 = JG+LDIN 


38 © J7+JDIM 


CALL JACK(Q,Q692) Q(52) A654) CIS) ,QCI6),QG57), 
Q(38) ,Q.0(34),Q,Q(51) ,Q(52) 


(Pace the array sterting 
locations to progres 
Jack) 


v.2 SUBROUTINE JACK (UL, Vi. ,PL, WC, DELRT DAPI ,DELY) ,DYP) ,CONTP ,WOQ,01.Q,VLQ,PLQ) 


This is the asin caleulation routine. Subroutine JACK inpute the eajority of the pregrem contre! end aede! 
peremeteze. The pregren tapute che terrein vie cubroutine MOUTMR, calculates and cote initial eenditions end then 
entere the time leep thee ealovletes the vector components of the vind speed and the eurface leyer heights. At the ond 
of each time leop, the program checks to detersine if printing or «pe output of che vind field ts desired for the 
present time step. When the tine loop hee been completed, the program loops to che next problem in sequence if present. 


hak oh fs 


Jun 
UL(L,J,K) 


VL(E,3 4K) 
PL(1,J,R) 
ULT 


ISKIP 
IUNIT 


PROGRAM VARIABLES. 


Array containing the inveree of (X(1+1) - X(I-1)] to minimise divisions in the tine loop. 
Arrey contetoing the inverse of (r(s+1) - r-1)) to minimize divisions in the time loop. 
Array contaiaing che inverse of [X(1+1) - X(i)]} to minimize divisions in the time loop. 
Array containing the inverse of (Y¥(J+1) - ¥(5)] to ainimize divisions in the cime loop. 
Time increment for the time loop. 

Totel time in seconds of the time loop. 

Array containing the terrain heights. 

Nuader of grid pointe in the y-axis. 

Nuaber of grid pointe in the x-axis. 

LU-1 

JSL-L 


ULQ(I+(J-1)*IDIMe (K-1)®IDIM®JDIM) - The u component of the wind speed times che layer depth where K = 2 tu 
the present time step and K © 1 is the past time step. 


VLQ(14+(J-1) #1 DLM® (K-1)®IDIN®JDIM) - The v component of the wind speed times the layer depth. 
PLQ(1+(J-1)*IDIMNG (K-1)*IDIM®@JDIM) - The surface layer depth. 

The initial u component of the wind speed or the mean wind speed depending on ISKIP(4). 

The initial v component of the vied speed or the direction depending on ISKIP(4). 

The initial height of che surface layer. 

Array containing the x-axis of the reference coordinate sysctea. 

Atrey containing the y-axis of the reference coordinate system. 


Array containing che time in minutes et which the wind field and layer heights are to be output victhin the 
time loop. Values ere in ascending order and the waximum is used as the stop value. 


Index of the first point in the x array at which uniform grid spacing occurs slong the x-axie. This value 
and IND, JST, and JND below sre output to tape for use in plotting the area of uniform grid spacing. 


Index of the first point in the y array at which uniform grid spacing occurs along the y axis. 
The ending indices in che x and y arrays at which uniform grid spacing stops. 
Number of iterative Cime steps between the recalculation of the time step increment OT. 


Stability factor for calculating the time step DT and maincain DT at a critical value. This value should 
be aa close to | as possible and still maintain program stability, Values of 0.90 to 0.95 generally 
saintein stability. 


Number of time steps between the spplicatione of s nine point soothing function. 


Reduced gravity factor equal co g(l-S) where g ts the acceleration of gravity (9.8 nfo’) and S is the ratio 
of the potential temperature at the cop of the layer over che potential temperature at the bottom of the 
layer. 


Program control options. Refer to the user instructions or the program listing for detaile. 


The Fortran logical cape unit for program tape output. If more than one reel is written, a | is added to 
IUNIT and @ reel te assumed mounted on the next unit. A second reel will be required only if many time 
stepe or many problems are being oucpuc. Refer to the user instructions or the program listing for the 
Cape format. 


The siniaum allowable layer depth. 

Array containing the x axis of the acandard WSMR grid in UTM (kilometers) coordinates. 
Array containing the y axis of the standard WSMR grid in UT (kilometers) coordinates.. 
Conecant used in the nine point smoothing function. 

Wueber of groups of UL, Vi end PL output to tape in any one problem run. 

Arcey containing record one of each file output to tape. 

Time of the last time step to process. 

Flag used to check if problems has tape output. 

Counter for the PRINT erray. 


Program etebility fleg where if eet to greater than sere the problem ia stopped and the program goes to 
the next case. 


Time step cumber at which tape output occere. 
Time at which cape output eccurs. 


Other program variables are weed for temporary eterage and tadexing. 


eee Re: 


SQ21 = 0. 70720678 (Program covetents) 
3 0.5 Bt Vis oe 


y . 


DM = 2 (Dimensional parametera used for 
IPIDIN © IDINCIDIN indexing through the arrays ULQ, 
LIDIN © IDIMN®JDIN VLQ, PLQ and 8GQ) 
TIEDIN = LIDINSKDIN ee 


mete pngran eae 


READ QUST! 
Semeliat containing LL, JL,X,¥, ISKIP ALK, PRINT ,ULT,VLT, 
PLT ,DTLI@A Cl , ISMPTH CHT, LUNIT 


(Convert from wind epeed end 
direction to vector u and v 


components) 


VLT © VLT#180.0 
ULB © ULT*SIN(VLTORAD) 
VLE © ULTOCOS(VLTeRAD) 
ULT © ULB 
VLE © VLB 


(Provide defoult values for input parameters not set.) 


fem ie 


NOUT(1) © NBUT(1)+1 
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(Provide default values for rhe 2 and y arrays if not set.) 


X(I) © XPS(E)41.0n3 
sth 5 pai (Convert te meters) 


(Determine coordinate indices of the area of waifere grid spacing vithin the grid.) 


CALL MISC(JST, 0D,¥, JL, JL) 


(laput che terrain heights.) 


eb. 


JL,2,%) 
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v 
(laitialise output tape 1f requested.) 


MBUT(2) © [sKIPC3) 
wiuT(3) = [DIM 
WUT (4) © IDIN 
1TaG) © 1 


{Output record 1 of this file) 


i usiglga 


aa 


i oe 


GIRBEEES © 1b O/ (aided )-0tt-12) 
GEPTEL) © 0 Of lattetd-oftn) 


DS 200 Jo2, 318 


DELYI(I) © 2.07 (¥(341)-¥(3-1)) 
DYPI(J) © 2.0/(¥(I+1)-¥(5)) 


(Output the initial iapus peremecets fres the Mewelist QLst!.) 


(Qatput che initial wind fielé, layer height and cerrein dete.) 


(Taitielize constants for time loop.) 


LIDINP = 1-1BL2 
WIDIMM ©@ -1-18L2 

Ist = 0 
Neo 


(Start time loop.) 


(Time etep counter) 


(Recalculate DT, time etep inccenent.) 


Pee = 1.0/PLQ(1) 
CMR © SQRT((VLQ(1)*PKR) 9% 24 (VLQ C3) @PKRRDO*2) CSQRE (GI MPLQ(1)) 


DT © UTLMDA*XMIN/XMAX 


oo ama 


(accumulate ¢ise.) 


(Solution unstable) 


(Evaluate UL, VL end PL ot each grid coordinate for thie tine otep.) 
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de teint enemies 


Caleulete relative indices fer the 


etere parameters that ere constant 
for the J dinension. 


JPL2 © J*IDINGIL+IBL2 
JM12 © JPL2-IPIDIN 
312 © JPL2-IDIN 
JZ = LIDUNEIL2 


DELY © DELYJ(J) 
DYP © DYP3(J) 
DYN © DYPJ(J-1) 


JPL2 © JPL2+1 
MLZ © M21 
312 @ Ji2el 
Ji3 @ JLW 


UL(I,3,1) 10 equivalent to 
ULQ(J12) and 
UL(I,J,2) © ULQ(J13) 
UL(I42,J,2) = ULQ(J1241) 
UL(T,J%1,)) = ULQ(IPI2) 


(The following equations do not reflect the form of the model equations 
due to slterations for calculation efficiency. Parte of the equations 
have been separated to eliminate redundant calculations and divistons, 
multiplications and array indexing have deen aiutal zed.) 


ULB © ULQ(SI2) 
VLB © VLQ(J12) 
GLAP © ULQ(J12¢1) 
VLXP © VLQ(J1241) 
PLKP © PLQ(J12¢1) 
ULEM ©. ULQ(Ji2~-1) 
VLEM © VLQ(312-1) 
PLUM © PLO(SI2Z~-1) 
ULYP = ULQ(SPI2) 
VLYP © VLQ(JP22) 
PLYP = PLQ(JP12) 
ULYN = ULQ(JH12) 
VLYM © VLQ(JH12) 
PLY © PLQ(JM12) 
ULXYP © ULQ(sP12¢1) 
VLEYP © VLQ(SPI2+i) 
PLXYP © PLQ(JPI 202) 
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ULXYN © ULQ(JMI2-1) 
VLXYM © VLQ(JM12-1) 
PLXYM © PLQ(IML2-1) 
VLONP © ULQ(JPL2-1) 
VLINYP © VLQ(JP32-1) 
PLXNYP = PLQ(JPI2-1) 
ULXPYM © ULQ(JMI2¢1) 
VLXPYM © VLQ(JM2 201) 
PLXPYM = PLQCINI 2¢2) 
PLBL = 1.0/PLB 


PLXPI © 1.0/PLXP 
PLXML © 1.0/PLE 
PLYPL = 1.0/PLYP 
PLYML = 1.0/PLYN 
PLXYPL = 1.0/PLKYP 
PLXYMI © 3.0/PLXYM 
PXNYPL © 1.0/PLEMYP 
PXPYMI = 1.0/PLXPYM 


AID = WCQ(J12-1BL2) 
AG3 © HOQ(IP12-16L2) 
A79 © HCQ(JM12-1BL2) 
A2k © HGQ(JL2+1DIMP) 
AIS © HCQ(IPI2Z+1 IDIMP) 
A76 © HCQ(IM) 241 JDIMP) 
A77 © HGQ(IP1 21 DIM) 
ATE = HGQ(JLZ+1 DIMM) 
ABO © HGQCIMI 241 DIMM) 
Al = 2. 0%DELKI(T) 

A2 © 2.0*DELY 

AlO = ULXP*PLXPL 

Al2 © ULXP*AlO 

ALL © ULB*PLBL 

Al4G @ ULDSALL 

AS = PLXPSPLXP 

Al6 © PLBSPLB 

A56 © 0, SDELY 

A46 © ULXYP*PLXYPI 

AL? © AGO*VLXYP 

ASE © ULXPYM*PXPYMI 
ALG @ ASL@VLXPYM 

A4D © ULYP4PLYPL 

Ald © ASS*VLYP 

AS3 © ULYN*PLYNE 

A20 = ASI*VLYM 

a6 @ A19-A20 

A33 = (PLXP+PLB)OG2 
A22 = A21-A39 


(Part of second-order terme common to UL, VL and PL) 
¥P = OXPL(2)4( (Al 2-A1l49G2%(A3-A1 6) )+AI3*%A22) A560" (AL 7-AlBtA6) 
» ff) Equation (4-20) in Section 4. , 


AG © ULXM*PLANI 
A28 © ULXM*AS 

A29 © PLXM*PLEN 
AGT © ULXMYPOPXMVPE 
A3O © AG7TEVLXMYP 
AS2 © ULXYM*PLXYMI 
AS) © AS2°VLAYN 
A‘O © (PLAPLXM)#C2 


(Part of second-order terms common to UL, VL and PL) 


9 © OXPE (1-2) 4¢ (41 4-A284628 (Al 6-A29) )4A40* (A39-A78) 4A56% 
(A6*A30-A32) 


© (f°) Equation (4-20) in Section 4. 


AS6 © 0, S*DELKI(I) 
A34 © ALOCVLEP 

A3S © ABSVLEX 
AIS © ADO-A35 

A? = VLYPOPLYPI 
A36 © VLYPCA? 

AGS © VLEOPLat 

A}? = ASSOVL 
A}@ © PLYPOPLYP 
AGS © C20 (PLYPOPLS) 


(Part of second-order terue common to UL, Vi and PL) 


GP © AS49(A17-03900A25) 4DVP* ( (436-0374 20 (038-016) 4A 4° 
(a43-A39)) 


¢ 
© (8) Equation (4-21) tn Section 4. 
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AS0 © VLDMFLINIL 
ab] © VL@OASO 

a42 © PLYNOPLYN 
A23 = G24(PLOPLID 


(Wert of eeccad-order term common to UL, YL end PL) 


GM © A54%(A25A16-A31 )ODYM* ( (A37-A6 14624 (Al 6-442) 44294 
{4}9-A79)) © (g-) Equation (4-21) in Section 4+. 


ASS © ASOPASS 
aS @ Al2-A28 
430 = Aj~A2S 
AS? = ASGPASO 
A26 © GlOPLB 
A24 © A2L-A16 


(Birat-order term of u component Equation (4-)5)) 
BL © OT*(A55*(ASHC22ASB) AST*AGTALOPAIFA 26) 


ATL © VUXYP-VLEIYs 
ASS = VLYP-VLYM 
AS = ALG*PLAl 
ALS © (ULXP-ULS)*OXPI(1)+a56* 
(A1144S9) 
(Part of second-order tera of u component, Equation (4-17)) 
BIP = (ALOFALL) @FP+AL3®(A33-0. S*(Al 2°PLXPI+A6)) 
AT2Z © VEXMYP-VLAYN 
ASS © (ULB-ULXM) *DXPI(I-1)¢ 
A56*(AS9%A72) 
(Pact of eecond-order term of u component, Equation (4-17)) 
B21 © (ALIL+A8) MFM+AKK* (040-0. 5* (AG+ABCA8) ) 
A27 © PLXYPOPLXYP 
AO2 © PLEMYP*PLINYP 
a4S @ ALLOVLE 
AJ} © ULXYP-ULINYP 
Ab? © ULXP-ULM 
ALS © (VLYP-VLB) *DYPCAS4* (47 AG?) 
(Pert of second-order term of u component, Equation (4-17)) 


BIP | 0.58( (AI+ASS) * (AS4* (AGO*ULAYP-AG7#ULXHYP+A 54628 (A27- 
AG2FAS8) ) +DYP* (Al 9-048) 4AGS4AS44 (AT 5-A7 7426) 4 (AGI 
AS1)@GP-(Al 9*PLYPICAGS*PLBI) “Al 5) 


A61 © PLXPYNOPLEPYH 
A63 = PLXYMOPLAYN 
A64 © ALB*PLII 
AT4& © ULXPYM-ULXYN 
A32 © DYM* (VLE~VLYN)+A5S6@ 
(A674A74) 
(Pare of second-order term of u component, Equation (4-17)) 
BE @ 0.59( (A456A50) #(A54% (ASPAS) MULKPYM-AS2°ULXYMEG 28 (ASB+ 
Ab1-A63) )+DYN® (A66~426)4A239A54% (A24+A7 6-A80) )+( Al 14453) 
OGH- (AG47A20*PLYMI) %A32) 
A770 = A36-A4l 
AGB © A36-A4Z 
A69 © AG3-A79 
(Taret-order term of v component, Equation (4-15)) 
Ch © DE®(AS5S#A2S0AS7*(A700C2%A68) A268A2*A69 
(Pace of second-order ters of v component, Equation (4-18)) 


CIP | 0.5¢( (VLEPOPLEPICAGS) OPPO (Al OFAIL) © (DKPI (1) ®(A34-A48)+ 
AS6* (VLAYPOVLAYP*PLXYP 1] -VLEPTMOVLXPYN® PKPYN[ +A 70+G 28 
ae )4AD3*AS6* (A75~A7 GFA69))~ (ADS*PLEPIFANG) 
Sal} 


(Pert of secend-erder tern of v component, Equation (4-18)) 


CAA © 0.50 ( (ASSOVLINOPLINE) OFM+ (Al 1408) © (DMP TE (T-1)*(A66~A8e 
VUIN)+AS6® (A7OOVL IY OVLINYPOPENYP I -VLLYROVLE YM¢ PLE YNI¢ 
saa 463) )PR40°AS6* (4699877 ABO) )~(AG4tA)5*P LMI) 
ears 

B66 @ A65*065 


(Peet of second-order tera ef v component, Equation (4-18)) 
CIP © (A7OR43) CPOA1S9(A63-0. S*(A7%A74066)) 
(Pert of second order tera of v component, Equation (4-16)) 
CH © (A450AS0) 0003 2¢ (025-0. 5° (A660K 900050) ) 

AGO © ASS*AG7ORS70059 
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te velocity conponest , Equation (4-14) 
WLQEIII) © VLQ(I12)~B1tA94 (AL* (27-821) 4A2* (BIP-BIM)wC) CDELLI (1) 
*426*A60) 


(wv velocity component (a, ,)) Equatson (4-15) 
VUQ(QEI) © VLQ(I12)-C14A9(A1* (C2P-CAN)4A2* (CIP-C 3M) Gh MDELY*AGS® 
A60) 


(ayer height (9, .)) Equation (4-16) 
PLQ(IL3) © PLQ(312)-DTAAGOPAD® (Al ® (FR-FM)+A2Z* (GP-GM) ) 


{ YES 


— oe 


(Set up the boundary conditions by aetting the ourermost roves 
end colyans equa] to the second to qutermost rows and columns) 


UL(2,J,KDIM) © UL(2,J,KDIM) 
VL(t,J ,.KDIM) © VL(2,J,RDIN) 
PL(L,J,XDIN) « PL(2,3,KRDIN) 
ULCLE,J,KDIM) = UL(LUN,J,KDIM) 
VU(LL J ,KDIN) © VLCLLM,J sia) 
PL(LL,J,KDIN) © PL(LUM,J,KDIM) 


UL(1,1,KDIM) © UL(1,2,KDIM) 
VLCE, 2 RDIM) © VL(2,2, RDI) 
PL(I, 1, ROIN) © PL(L, 2, RDU) 
ULC ,JL,ADIM) = UL(1, JUN, KDIN) 
VUCT, JL, RDU) © VEL JUN, KDI) 
PLCL,JL,KDIN) © PE(1, JUN, KDI) 


PLQ(L) = PLQ(I+1J0IM) 
Is 
PLQ(L) > 20000 >>2E2 
? 
: aaa 


UL(L, J, RDiMeL) © ULCT J KDIM) OB 10 (UE (161, J KDIM)CUL(I-1 J KDIN}S 
UL(1, 348 RDIMDUL (I Jol KDIN)=4 0%. (I,J .KBIN) ¢Ci *(UE(t-2 J-2, 
DIN) FUL (2-1 361 , RDIN) OUT. (Tol, J-1 , KDIN) SUL (141,341 .KDIM) 6.08 
UL(T,3,6DIN)) 
VE(E J, KOIN-3) © VLCR 2, SKIN) OB) (VL (14) ,J , RBIM)OVL (1-1, 5, RDIN)S 
VLC Jos KDIN) OVE (E, 2-2 DIM) -4. OOVL (1,3, KBIN) )4C1IO(VL(I-1,3-1, 
RDIND SUL (Lo) Jo) ADIN) SVL (11 ,J-1 , MDIW)OFL (161 , 341 , RDIN)-6.00 
VL(2, J, XDIN) 


(Posh peremetera down tins steck Sates encething te called for.) 
YES 

4 

[v6 00 t=t, 13000 . 

(Doves tarrein penetrate leyert) 
ULQ(L) & ULQ(I¢I JDIM) 
VLQ(Z) © VEQCI +I DIM) 
(ie problem unstable?) 
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PL(I,3 ,KDIN)) 


PL(1, 3, RDIN-1) eABLK 


PL(I,3,KDIM-1) < ABLK 
? 


is 
PL(I,J3,KDIM-1) > 20000 
? 


ves LELAGeN 


no 


620 CONTINUE 
asia 


- 


NO 


eas ea eet Bere | 


a ae 


(Apply boundary condittons co emoothed values.) 


+ tees 


Dd 640 Jel JL 


ULCL,J,KDIM-1) = UL(2,J,KDIM-1) 
VL(L,J,KDIM-1) = VL(2,J5,KDIM-1) 
PL(1,J,KDIM-1) © PL(2,J,KDIM-1) 
UL(LL,J,KDIM-t) © UL(LLM,J,KDIM-1) 
VLCLL,J,KDIM-1) © VL(LUM,J ,KDIM-1) 
PL(LL,J,KOIM-1) © PLCLUM. J, KDIM-1) 


wot 


C1.(T,1,RDIMN-1) © UL(1,2,KDIN~1) 
VL(E, 8, KDIM-L) © VL(S,2,RDIM~1) 
PL(L,1,KDIM-L) = PL(1,2,KDIN~1) 
UL(L,JL,KDIM-1) © UL(T,JLM,KRDIM-1) 
VUCE, JL, KDIN~L) = VL(E,JUM, KDIN-1) 
PL(1,JL,KDIM-1) = PL(E,JUM,KDIM-1) 


: PL(I,3,KDIM-1) = PL(1,J,KDIM6)+B2¢(PL(I¢1 J, KDIM)+PL (1-1, ,KDIM)+ 
7 PL(L,J4+1 , RDEM)OPL (1, J-2 ,KDIM)-%.O*PL (1,5, KDIM))¢C14(PL(I-1,J~1, 
KDIM)¢PL(1-1, 541 ,RDIM)@PL( 14] ,J-1, KRDIM)+PL (11, Je) ,KDIM)-4.08 
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(Check for unstable sulut ion) 
YES 


(Output wind field and leyer heights if this te « selected time.) 


1s 
CIM < PRINT (IPR)*60.0 
? 


(Output wind field and layer heights.) 


WRITE N 
Meseage thet solution 
ie unstable 


Ke rsPs(i) 
TIME © TIM®.2777777E-3 
LINES © 60 


PLB © 1.0/PL(2,3,3) UL(1,3,2) © ULB 

ULB @ W(2,J,1)°PLE VL(1,J,2) © VLB 

VLE © VL(1,J,1)¢PLE PL(I,3,2) © PLS 

PLE @ PL(T,J,1)00G(1,9) LINES © LIWES+1 
voat = 0.0 


YES 


VUXP @ VL(it} J,3)/PLOIe1, 1,1) 
VEXM © VL(I~3,5,3)/PL(1-1,3,1) 
ULYP © ULC(E,J+l,3)/PL(T,J45 i) 
ULYM = UL(L,J-1,2)/PL(t,3-3,1) 
VORT = (VLXP~VLXM)*DELRI(1)~ 

CULXP-UL IM) 4DELY ICI) 


yes 


WRITE New heading, 
new page 


(Convert from vector components to 
speed and direction) 


YES 
CALL DYDIA(DLB, VLB) 


I 


ves 
= 
ue 
(Qucput wind. field gad layer heights co eave tape.) 


ISAVE(JPR-1) = ¥ 
TSAVE(IPR-2) © T1M/60.0 
CALL OUTPT(IUMIT ,LJDIM, ULQ(IJDINe!) .NEGS) 


(Check fer lest time step) 


aie, ot eg eee seg eee ae! ee 


(Close output cape file) 


| 


(Write tape output file conatents) 


WRITE (ISAVE(L) ,TSAVE(E), 
fet, {RECG) 
time steps and tise 


WRITE (MP=UT(I),I=3,10) 
OR umber of time stepa, ISKIP(3), 
LDIM,JDIM,LL,JL,IST, 
IND, JST, IND 


(Te thie the last or only case) 


Feet one 


0.3 SUBROUTINE MPUTWR(HC,LL,JL,X,¥) 


« Subroutine MOUTHR reade the terrain deta and the indices of the terrain data. The indices IST, IND, JST and 
; ° IND of the etarting and ending pointe on the x and y axee are input firet. The terrain height data are then input. 


i‘ Aveas of the grid that do not contain terrain height dete and automatically the firet (lest) three rows and columns 
; are filled by the program. These arese are filled vith terrain heights thet slope down to the mininua height at the 


boundary edges. 
QEaD IST, IND, JST, IND 
Terraia hetght indices 


HG(1,J) © XMIN 
WG(L,J¢SL-3) © XNIN 


Leo ep 


4 re) 


XML © 1.0/(X (IML) -K(iM3)) 
IL © IMl-1 


i DP SO I=4, IL 


WG(L,J) © CHGCIML, 3)-XHIN)@XM2¢XNIM 


: a | 
4 | 
| 

| 

| 

| 

| 

! YES 


XM) © 3.07 (X(LL-2)-x (1N3)) 
IL © LL-3 
IS = Tel 


WG(1,3) © (WCC DNS, 3)-tIN) CANQeITN 


tana eieel 


< 
r 
wv“ 


= 


wed = 1.0/(¥C1M2)-7¢(3)) 
tL = im2-1 
IR © UL-3 


| 
| 
J 


| 
{ 
| 
| 
| 
| 
| 
| 
4 


WA © 8.0/(¥(SL—2)-¥ (14) 
tL © JL-3 
{iS © INéel 


WO(1,3) © (MGCL, U4) -WIROeZeRNIA 


D-24 


80 CONTIWUE 


cote oe ae eae aod 


D4 SUBROUTINE GUTPT (LUNIT ,WWRD ARRAY ,NOCO) 


Bubroutine GUTPUT vrites the wind field and layer height information cto save tape. This routine vees the 
UMIVAC 1108 MTRAM eubroutines for all output. 


CALL NTRAN(IUNIT, 2 ,NWRD, ARRAY, ISTS) 


(NWRD words of binery data are vritten to tape IUNIT from the area ARRAY. 
The status of the operation is stored in ISTS at the completion of the 
operation.) 


(Loop if operation is not complete cto prevent overlap processing) 


(Branch {f write vas successful) 


(Branch if end of tape is detected) 


pcm ARE oN Ae Semi a nt Eee mi te ee 


(dvanch if parity or ndrévere error occurred) 


(Weitte a meusege thet an uareceguisable error code hae occurred) 


(Bachepace oma record aad write an end of file vu vuspur tepe and 
factrement the tape unit designator IUNIT. The last record iy 
then revritrem on the second tape.) 


CALL MTRAM(IUNIT,22,7,-1,9,11) 


(Rackapece to etart of record and retry the write. [f unsuccess- 
ful after S tries, return to calling program write error 
condition set.) 5 


(Set unrecoverable error, W6CH=1) 


9009 = 
cc) 


2.5 SUBROUTINE MISC (TIST IID, XY, LLIL LLILN) 


Subroutine MISC determines the starting and ending indices on the x or y axis that determines the area withia 
the grid ever which vailorm grid epacing occurs. These peremetere LIST and [JM for both axes are output to the save 
tape for plotting of the grid dependent parameters. 


OO 50 el ,LLJUN 


ws al 


ity 


(Bading parameter occurs the first tine the Jiffcrcace S22 
greater than the siniaun.) 


(Hag the minimus and [JST been determined?) 


iat 
XIN © DIF 
13ST © 0 


(Is the sininum equal to the last difference?) 


YES 


(Set starting tadex to leat sintmum detected.) 


LIND@LLJL 


D.6 SUQROUTINE UVDIR(X,¥) 


; Subroutine UVDIR converte the vector componente u and v of the wind speed into the vind speed and direction. 


DIR © 270.0-ATAN2(¥ ,X)#57.29582 


DIR-VIR- 300.0 
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This report describes the development and documentation of a computer algorithm, 
contalning a two-dimensional shallow-fluid model, for calculating the wind field above 
comiex terrain. The optimum finite-differencing procedure employed.'for numerical 
solutions of the’algorithm 's a Lax-Wendroff scheme using the grid points and two time 
tevels in combination with a nine-point low-pass fiiter. 

A compretiensive computational program, using an Isolated symmetri¢al mountain, was 
carrled out to provide guidelines as to the nature of two-dimensional solutions of the 
shallow-fluld equations for the wide variety of Initial conditions encountered In the 
atmosphere. The results showed that the flow patterns could be divided into four 
major categories: Subcritical without hydraulic jumps; supercritica’l without upstream 
waves; critical with hydraulic jumps and wind-direction reversals; critical with 
hydraulic jumps but without wind-direction reversals. For the subcritical and super- 
critical flows, Initialization procedures do not appear to pose a problem. for the 
critical flows, care must be taken In the selection of Initialization procedures. 

Comparisons of calculated wind fleld patterns with recent detailed observations of 
J wind circulations above comilex terrain show excellent qualitative agreement in the 

limited cases avaliable for analysis. Additionally, the computer algorithm for the 
two-dimensional model, when applied to the terrain at White Sands Missile Range, gave 
results that were consistent with the limited observations ava'iable for two exanple 
situations. The computer program containing the two-dimensic | shallow-fluid model, 
written In Fortran V language and designed for use on a UNIVAC 1108 machine, is fully 
documented In the appendices to the report. 
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